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STRUCTURE-PRESERVING FINITE ELEMENT METHODS FOR STATIONARY MHD 

MODELS 


KAIBO HU AND JINCHAO XU 


Abstract. In this paper, we develop two classes of mixed finite element schemes for stationary magneto¬ 
hydrodynamics (MHD) models, one using the magnetic field B and the electric field E as the discretization 
variables while the other using B and the current density j as the discretization variables. We show that the 
Gauss’s law for magnetic field, namely V • B = 0, and the energy law for the entire system are exactly 
preserved in both of these finite element schemes. Based on some new basic estimates for H^{div), we 
show that both of these new finite element schemes are well-posed respectively under some conditions for 
the B-E formulation but no conditions for the B-j formulation. Furthermore, we show the existence of 
solutions to the nonlinear problems and the convergence of Picard iterations under some conditions. 


1. Introduction 

In this paper, we will develop structure-preserving finite element discretization for the following 
stationary incompressible magnetohydrodynamics (MHD) system: 

(i.ia) {u • V)ii — — Sj x B Vp = /, 

IhfO 

(i.ib) j - X B ^0, 

■^rn 

(i.ic) V X = 0, 

(i.id) V • B = 0, 

(i.ie) V ■ u — 0, 

where the Ohm’s law holds; 

( 1 . 2 ) j = E + u X B. 

Here u is the velocity of conducting fluids, p is the pressure, B is the magnetic field, E is the electric 
field and j is the volume current density. Dimensionless parameters Re, Rm and S are the Reynolds 
number of fluids, magnetic field and the coupling number respectively. 

In the study of magnetohydrodynamics (MHD) system, it is well-known that the Gauss’s law for the 
magnetic field, namely V • S = 0, is an important condition in numerical computation of MHD system 
[5, 8]. Nonzero divergence of B will introduce a parallel force, which breaks the energy law. In our 
previous work Hu, Ma and Xu [12], we proposed a class of structure-preserving and energy-stable finite 
element discretizations that exactly preserve the magnetic Gauss’s law on the discrete level for the time 
dependent MHD systems. The goal of this paper is to extend such discretizations to stationary cases. 
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Such a discretization is however not straightforward as the time-dependent and the stationary systems 
have different structures. In the time-dependent problem, the Faraday’s law reads; 

dB ^ ^ 

-f V X £; = 0. 
at 

In [12], we chose to keep the electric field E and use the iF(curl)-conforming finite element space for 
E and iF(div)-conforming finite element space for B to discretize the above Faraday’s law as follows: 


Af 


-f V X E" = 0. 


This implies that V • B^ = 0 holds for all n > 1 as long as it holds for n = 0. 
In the stationary case, the Faraday’s law reads: 


V X = 0. 


In this case, we can not directly apply the technique used in [12] for the evolutionary case to preserve 
the Gauss’s law V • B = 0 exactly on the discrete level. Instead we will treat the Gauss’s law as an 
independent equation in the whole MHD system and we will then introduce a Lagrange multiplier to 
appropriately enforce this law on both the continuous and the discrete level. 

The idea of the use of Lagrange multiplier itself is not new (see Schotzau [15] and the reference 
therein) and the novelty of our approach here lies in how this technique is used in combination with the 
techniques developed in [12]. In Schotzau [15], a magnetic multiplier r G is used to impose 

the Gauss’s law in the following way: 


B • Vs = 0, 


Vs e 


which does not guarantee that the Gauss’s law holds strongly (namely V • = 0 point-wise in the 

domain) in the corresponding discrete case. The main difference in our approach is that the Gauss’s law 
will indeed be preserved on the discrete level strongly by using appropriate finite element discretization 
of B so that Bh is iV(div)-conforming. The finite element de Rham sequence as studied in [1, 11, 4] 
plays an important role in the construction and analysis in our paper. 

MHD equations admit many different variational formulations which lead to different mathematical 
properties and numerical efficiency on the discrete level. In most existing literature, variables E and j 
are eliminated to reduce the size of the corresponding discretized problems. In [12], we demonstrated 
that it is advantageous to keep E and use it as an independent (or intermediate) discretization variable in 
appropriate finite element space. Indeed, this approach may lead to larger discretized systems, but these 
systems have better mathematical structure and may be solved, as illustrated in [12], more efficiently 
than the corresponding smaller systems derived from traditional schemes by eliminating both E and j. 

In this paper, we continue and extend this study for the stationary problem by keeping either E 
or j and discretizing them with appropriate finite element spaces. We will refer the corresponding 
discretizations as B-E and B-j formulations respectively. 

For simplicity of exposition, we will use the following homogeneous Dirichlet boundary conditions 
for B-E system: 


M = 0, 

B ■ n = 0, 
E X n = 0. 
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According to the Ohm’s law that j = E+ux B, the above boundary conditions are obviously equivalent 
to 

u — 0, 

B ■ n = 0, 
j X n — 0. 

These conditions turn out to be all essential boundary conditions. The extension to non-homogeneous 
essential boundary conditions or non-essential boundary conditions is straightforward and standard and 
the relevant details will not be given in this paper. 

The rest of the paper is organized as follows. In §2, we present the notation and basic finite element 
spaces used in the discussion. §3 demonstrates basic estimates for H^(divO) functions, including regu¬ 
larity and the discrete Poincare’s inequality. §4 extends the B-E formulation to stationary problems. In 
§5, a new formulation based on B and j is studied. Concluding remarks and the comparison of B-E 
and B-j formulations are given in §6. 


2. Notation and basic finite element spaces 


In this section, we will introduce some basic Sobolev spaces and their corresponding finite element 
discretizations that will be used in the rest of the paper. 

We will assume 17 to be a bounded domain with Lipschitz boundary. For the ease of presentation, 
we further assume 17 is simply connected (contractable), so that there is no nontrivial harmonic form. 
Similar discussions as in [1] hold for more general cases. Using the standard notation for inner product 
and norm of the space 


/ u ■ vdx, ||u||:= ( 

[ Iwpdx^ 

/n ' 

\Jn J 


1/2 


we define the following space with a given linear operator D: 


:= {u e L ^( fl),Dv S L^(17)}, 


and 

Hq{D, 17) := {u S H{D, 17), tov = 0 on 917}, 
where tu is the trace operator: 


tov := 


V, D = grad, 

V X n, D = curl, 

V ■ n, D = div. 


Here iJ(grad, 17) is a scalar function space, while iT(curl, 17) and iJ(div, 17) are for vector valued 
functions. And we often use the following notation: 

Lliil) := G L2(17) : j v = 0 


When D = grad, we often use the notation: 

iT^(17) := iF(grad, 17), iTo(17) := iFo(grad, 17). 

For clarity, the corresponding norms in H{D, 17) are denoted by 

\\u\\l= ||«|p+||V«f, 

ll^llcurU= ||-Ff + |!VxF||2, 
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|C||^w:= 


.|C'ir + ||V-Cf. 

We will also use the space with norm || -Ijo.p given by HuHg p= For a general Banach space 

Y with a norm |j -Ijy, the dual space Y* is equipped with the dual norm defined as 

{h,y) 


|fi.||y:= sup 

O^yeY 


\y\\Y 


For the special case that Y 
which is defined as 


iFo(n), Y^* = H ^(n) and the corresponding norm is denoted by IHI-i, 


|/||_i:= sup 


if^v) 

llVttll 


We will use Ci to denote the constant in the following inequality, which is a consequence of Sobolev 
imbedding theorem and Poincare’s inequality: 

||u||o,6< CillVull, V'ueiFo(n). 


Since the fluid convection usually appears in the following discussions, we introduce the trilinear 
form 

L{w;u,v) := -[((u; • V)u,v) - {{w ■ V)r;,M)]. 

When lu is a known function, L{w, u, v) is a bilinear form of u and v. This will occur in the Picard 
iteration, where w is the velocity of the last iteration step. 

Let This, Si triangulation of and we assume that the mesh is regular and quasi-uniform, so that the 
inverse estimates hold [6]. The finite element de Rham sequence is an abstract framework to unify the 
above spaces and their discretizations, see e.g. Arnold, Falk, Winther [1, 2], Hiptmair [11], Bossavit [4] 
for more detailed discussions. Figure 1 shows the commuting diagrams we will use. Electric field E, 
magnetic field B and the multiplier r will be discretized in the last three spaces respectively. Figure 2 
shows the finite elements of the lowest order. 


iLo(grad) iTo(curl) Tfo(div) 

I nr' I nr |n 

iTo^(grad) 7Tg^(curl) i?o(div) 

Figure 1 . Continuous and discrete de Rham sequence 


0 

h 



Figure 2. DOF of finite element de Rham sequence of lowest order 

As we shall see, H (div) functions with vanishing divergence will play an important role in the study. 
So we define on the continuous level 


iLo(divO,0) := {C S iLo(div,0) : V • C = 0}, 
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and the finite element subspace 

i/^(divO, 17) :={Ch& -ffo^(div, 17) : V • = 0}. 

We use Vh to denote the finite element subspace of velocity Uh, and Qh for pressure ph- There are 
many existing stable pairs for Vh and Qh, for example, Taylor-Hood elements [10, 3]. T/q (div, 17) and 
Lq ^(17) are finite element spaces from the discrete de Rham sequence. For these spaces we use the 
explicit names for clarity, and use the notation Vh and Qh for the fluid part to indicate that they are 
usually different from Hq ^(17)^ and Lg ;j(17) in the de Rham sequence. 

There is a unified theory for the discrete de Rham sequence of arbitrary order [3, 1, 2]. In the case 
n = 3, the lowest order elements can be represented as: 

K ^ a V 2 AP A ViAP a VqA^ 0 

K ^ V2AP a ViAP a a VoA? 0 

K ^ V 2 AP A V^JA^ A ViAP a VqA^ 0 

K ^ ViAP a V^A^ a a VqA^ 0 

The correspondence between the language of differential forms and classical finite element methods is 
summarized in Table 1 . 

To link the finite element spaces, below we will require iTg (curl, 17), iJg (div, 17) and Lg ^(17) to be 
in the same sequence. 


k A* (17) Classical finite element space 


0 VrA°{T) 

1 VrA\T) 

2 VrA^{T) 

3 VrA^iT) 

0 V-A°{T) 

1 V-A\T) 

2 V-A^{T) 

3 V-A^{T) 


Lagrange elements of degree < r 
Nedelec 2nd-kind H (curl) elements of degree < r 
Nedelec 2nd-kind H (div) elements of degree < r 
discontinuous elements of degree < r 

Lagrange elements of degree < r 
Nedelec Ist-kind iT(curl) elements of order r — 1 
Nedelec Ist-kind 7T(div) elements of order r — 1 
discontinuous elements of degree < r — 1 


Table 1. Correspondences between finite element differential forms and the classical 
finite element spaces for n = 3 (from [1]) 


As we shall see, it is useful to group Vh x iTg (curl; 17) x iFg (div; 17), and define 

Xeh := Vh X iTo^(curl,17) x iTo'‘(div, 17), 
and group Qh x Lg ^(17) to define 

Yh-.= QhxLl^{n). 

For the B-j formulation, we group the spaces to define 

Xjh ■= Vh X iTQ(curl, 17) x iTQ(curl, 17) x iTg(div, 17). 

For the analysis, we also need to define a reduced space, where jh and cTh are eliminated: 


X,h ■■= Vh X iTo^(div,17). 
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In order to define appropriate norms for the B-j formulation, we introduce the discrete curl operator 
on the discrete level. For any Ch S Ffg (div, 17), define Wh x Ch € TFq (curl, 17) as: 

(V,, X Ch, Fh) = {Ch, V X Fh), yFh e iJo^curl, 17). 

And for any Wh G Ffg (curl, 17), we define Wh ■ Wh G Hq (grad, 17) as 

(V/j • Wh, Vh) = -{wh, Vu,i), yvh G iFo (grad, 17). 

We define P : L^(17) —>■ TJq (curl, 17) to be the projection 

{F<j,,Fh) = {(j>,Fh), ^Fh G if^(curl,17),(^ e L\n). 

We further define |j-||(i to be a modified norm of iFg (div, 17) by 

11^^11?:= l|C'„f + ||V-C^f + ||V;,xC;.f. 

And |j 'lie for iJg (curl, 17) is simply the norm: 

\\Fh\\l-.= \\Fhf. 

There are some motivations to define such a stronger norm for Hq (div, 17) and weaker norm for Hq (curl, 17) 
space. One technical reason is that we want the nonlinear term V x {uh x Bh) to be bounded in some 
proper discretization. But generally Uh x Bh may not belong to iFg (curl) for Uh G 7Tg(17)^ and 
Bh G iFo(div, 17). So we choose to move the curl operator to the iTg (div) test function in the varia¬ 
tional formulation to get (uh x Bh,'^h x Ch)- Therefore we add the weak curl norm to iJg (div, 17) 
space. Another motivation can be seen in the energy estimate: on the continuous level, the energy esti¬ 
mate contains j — x B, but not V x j. So it is natural to use norm for the discrete variable 

jh- 

Now we define the norms for the various product spaces. For Yh space, we define 

||(<7,r)||^:= Ikf+llrf. 

For the other product spaces, we need to use different product spaces and different norms in the two 
different discretization schemes. For the B-E formulation, we define 

\\{u,E,B)\\\^-.= ||«||?+||£;||Li+I|S||L. 

For the B-j formulation, we define 

\\{uh,jh,(Th,Bh)\\xy.= \\uh\\l + \\jh\\l+\\o-h\\l+\\Bh\\l, {uh,jh,(Th,Bh) G Xjh, 

and 

\\iuh,Bh)\\l:.= \\uh\\j+\\Bh\\l {uh,Bh) G X,h. 

■^3 

3. Estimates for divergence-free vector fields 

In this section, we will establish some new regularity results for the strong divergence-free space 
iTg (divO, 17) which will be used for our forthcoming analysis. The main ingredients used in our analysis 
include some regularity results for the space Z := iF(curl, 17) n iFo(divO, 17) (c.f. [11, 15]), and for the 
space 

XI ■= {tu e iFg (curl, n) :Vh-Wh = 0} 

(c.f. [11, 15]), together with some appropriately defined “Hodge mapping” {Hd below) that connects 
iTg (divO, 17) with Z. 

We first give a preliminary result based on Hodge decomposition: 
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Lemma 1. 

V X Z = i7(div0,17) = V X H{cur\, 17). 

Proof. From the Hodge decomposition for L^(17)^; 

= ViFi(17) + V X iFo(curl,17) = iJ(curlO, 17) + iFo(divO, 17). 

Here iJ(curlO, 17) := {F e iF(curl, 17) : V x F = 0}. 

Therefore 

iF(curl, 17) = n iT(curl, 17) 

= H (curio, 17) n iF(curl, 17) + Hq (divO, 17) n iJ(curl, 17) 
( 3 . 1 ) = iF(curlO, 17) + Z. 

This implies 

iF(divO, 17) = V X iF(curl, 17) = V x Z. 


□ 


We now define the “Hodge mapping” for TJq (divO) functions. Let Hd : TFq (divO) Z be defined 
by 

( 3 . 2 ) {\/X {HdBh),\7 X v) = {\/h X Bh,\7 X v), \lveZ. 

Due to the Poincare’s inequality of Z, || 2 :||< ||V x z\\ holds for any z G Z. Therefore ( 3 . 2 ) uniquely 
defines HdBh- 

From Lemma 1, we have V x Z = iF(divO). Therefore 

( 3 . 3 ) (y X {HdBh),w) = {\7h Bh,w) , Vm e iF(divO). 

In particular, choosing w — W x (HdBh), we see 

||V X iHdBh)\\< llVi, X Bh\\. 

In the following, we will use B to denote the continuous lifting of B^: 

B := HdBh. 

And He : > Z is the Hodge mapping for iTg (curl, 17) [11, 15], defined by 

V X {HeFh) = VxFh, ^Fh e XI 
We also use the notation F to denote H^Fh when Fh G X^. 

Lemma 2 (Approximation of Hd). If ft is a bounded Lipschitz domain in 

\\Bh-HdBh\\<h--+^\yh X Bh\\, 

for all Bh G iTg (divO, 17) with 0 < <5(17) < i. 

Proof. We define to be the bounded cochain projection to Hq (div, 17) [9] . Note that V • (^B^ — 

0 due to the commuting diagram. Therefore there exists <pi^ G X^ and the corresponding lifting 
0 := Hc4>h G Z such that Bh — — V x cfh — x ^ and 

(3-4) Uh - X cf>h\\= - nSi.BII, 

where the first inequality is from the approximation property of H^. 

From ( 3 . 3 ), we have 

(Vft X BhA) = y X By) = {By x <^), 
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and 


Namely, 


Thus 


{Bh, V X cl)h) = {Vh X Bh, (fJh) = (V/, X Bh, 0/1 - 0) + (B, V x ^). 

{Bh — B,Bh — = (V/i X Bh, c()h — 0). 

IIB;, - Bf = {Bh -B,Bh- n^B) + {Bh - B, - B) 

= (V,, X Bh, ct>h-^) + {Bh - B, n^B - B). 


By ( 3 . 4 ) and the interpolation error estimates 

||B - n5i,B|i< /i5+^||B|| i+,< X B||< hi+^W^h X BhW, 

we obtain 


{Vh X Bh,4>h - <A) 


and hence 

\\Bh- 


This completes the proof. 


</l5+^|lB;,-nLB||||V^xB;,|| 

< h-^+^ (\\Bh - B\\ + \\B - n^BIl) IIV,, X BhW 

< h^+^WBh - SIIIIV;, X B^IIW+^'^IIV,, X BhW^ 

< ^WBh - BfX BhW^W+^^W^h X BhW^ 

B||2< ||B-n;ji,Bf+/ii+2^||v^ X B^f. 


□ 


For nonlinear problems and their linearizations, it is technical to prove the boundedness of variational 
forms, and this often requires careful estimates of regularity. The nonlinear terms in the B-j variational 
form of MHD equations have the form {uh x Bh,jh), where Uh € Vh C Bq Bh S Bq (divO, fl) 
and jh e B^(curl, fl). 

Lemma 3. For Uh G Vh and Bh G Bq (divO, fl), we have the following bound: 


W^h X BhW< |jM/i||i||V/j X BhW- 


Proof From Lemma 2, we have 

||B;,-B|1<LH^||V;, XB;,||, 

where 0 < S < | is a positive constant depending on the domain. 
Then 


Wuh X BhW< Wuh X {Bh - B)W + Wuh X B||. 

For the first term, 

Wuh X {Bh - B)|| < lltt/illo.ooll-Bft - -B|l 

< h~^WuhWo,6-h^W'^h X BhW 

< ||m,,||i||V,, xB,,||, 


where the second inequality comes from the inverse estimates and the approximation results. 
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And due to the regularity of Z [11], we have 

\\Uh X B\\ < ||m/i||o,6|!-B||o,3 

— l|tt/i||i|lV/j X BfiW- 


This implies 


\\uh X Bh\\< ||M;i||i||Vft X Bh\\. 


□ 


Below we will use a positive constant C 2 to denote the bound; 

\\uhxBh\\<C2\\Vuh\\\\VhxBh\\, 

and therefore 

{Uh X BhJh) < C2||VM;i|||jV?, X BhWWjhWo- 

In the discussions below, we will need discrete Poincare’s inequality for I/q (divO, fl) functions. We 
note that the two dimensional case is given in [7], and the proof can be modified to adapt to the three 
dimensional case. We include a different proof here. 

Lemma 4. For Bh G T/q (divO, 17), we have the following discrete Poincare’s inequality: 

\\Bh\\< llVi, X Bh\\. 

Proof. Because V • Bh — 0, we can choose Eh G Hq (curl) such that 

V X Eh = Bh, and Vh ■ Eh = 0. 


From the discrete Poincare inequality for in [2], 

(3-5) ||L;i,||curi<||Vx£;^||=||B;,||. 

And we have 

{fJh X Bh,Fh) 


X BhW = sup 

(curl) 


(3-6) 


= sup 

FfeGffo^(curl) 


(B/,, V X Eh) 


Therefore combining ( 3 . 6 ) and ( 3 . 5 ), we get 


and 


IV7 D ll\ X Eh 

|V.xB,||>.-^-^^^ll. 

I!v„ xBJ|> IIBJI. 


□ 


4 . B-E FORMULATION 

The variational formulation and the corresponding finite element discretization based on B-E are 
similar to our previous work [12]. Since our focus is the finite element discretization, in this section 
we give detailed discussions on the discrete level. The same conclusion also holds for the continuous 
formulation with minor modification. 
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4.1. Variational formulation and mixed finite element methods. We propose the following varia¬ 
tional formulation for (i.ia)-(i.ie): 

Problem 1. Find (u, E, B) £ X eh i^nd (p, r) G Yh, such that for any (v, F, C) € X Eh ond (q, s) € 
Yh, 

L{u\ u, v) + Vn) 

Re 

( 4 . 1 ) - S{j X B,v) - {p,W ■ v) = {f,v), 

( 4 . 2 ) 5(j,F)-^(B,VxF)=0, 

■^rn 

( 4 . 3 ) ^(Vx£;,C) + (r,V-C) = 0, 

( 4 . 4 ) -(V-M,g) = 0, 

( 4 . 5 ) (V-B,s)=0, 

where j is given by Ohm’s law: j = E + u x B. 

We can verify various basic properties of the variational formulation: 

Theorem 1. Any possible solution of Problem 1 satisfies 

(1) Gauss’s law of magnetic field in the strong sense: 

V • S = 0. 

(2) the Lagrange multiplier r = 0 and ( 4 . 3 ) is reduced to 

V X £; = 0, 


which holds strongly. 

(3) the energy estimates: 

and 

Proof The Gauss’s law directly follows from ( 4 . 5 ). 

Now taking C = V x in ( 4 . 3 ), we get V x E = 0 and hence 

(r, V • C) = 0, VC G (div, 12). 

This implies that r = 0 since Lg hi^) ~ ^ ' ^0 (both on the continuous and discrete level). 

The first energy estimate follows by taking v = u, F = E,C = B in ( 4 .i)-( 4 . 3 ), and then adding 
them together. The second energy estimate follows from the following inequality: 

(/,«) < ||/||-l||V«||< ^||Vw|P + ii?e||/||^i, 


This completes the proof. 


□ 
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4.2. Picard iteration. We give the following Picard iteration schemes for solving nonlinear system 
(4-i)-(4-5)- 

Algorithm 1. Given find G Xeh and G Yh, such that for any 

{v, F, C) G Xeh and {q, s) G Yt, 



-f Vrt) 

Fte 

(4-6) 

-^(j:_ixB"-i,r;)-(p",V 

(4-7) 

^(jn-i,-F’)-7^(S",VxF)=0, 

(4-8) 

^(V X F”,C)-f (r",V-C) =0, 

^rn 

(4-9) 

-(V-«”,g)=0, 

( 4 . 10 ) 

(V-B”,s) = 0, 

where 

/i defined as: -f m" x B"“^. 


The divergence-free property, consistency and energy estimates can be obtained similarly. 

Theorem 2. Any possible solution of Algorithm 1 satisfies 

(1) Gauss’s law of magnetic field in the strong sense: 

V • = 0. 

(2) the Lagrange multiplier r" = 0, hence (4.8) reduces to 

V X £;” = 0. 

(3) the energy estimates: 

^||V«"f+5|lj^if= (/,o, 

2le 

and 

(4.11) 

Next we show that the Picard iteration is well-posed under some assumptions. For simplicity, we 
denote ^ = {u,E,B), rj = {v,F,C) and x = (p, r), y = {q,s). u~ and B~ are some known 
functions. 

Define 

aE{i,v) ■= 2 -(Vm, Vt!) 

IXq 

-\- S{u X B~,v X B~) + S{v X B~,E) + S{u x B^,F) 

+ S{E, F) - V X F) + ^(V X F;, C), 

and 

y) ■■= -(V • M, g) -f (V • B, s). 

Algorithm 1 can be recast into the following mixed formulation with (it, E,B,p) = (it", F", B",p"), 
u- = and B^ = B""!; 
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Problem 2. Given 6 G and i/i G Y^,find {u, E, B,p, r) G Xeh x such that 
( 4 . 12 ) aE{i,'n)+ b{ri,x) = {e,rf), ^ "H & Xeh-, 

( 4 - 13 ) b{i,y) = {i>,y), \fyeYh. 

Theorem 3. Problem 2 is well-posed under the assumption that B~ G and that 

( 4 - 14 ) Re < -^C^^WfWzl 

It goes without saying that the assumption ( 4 . 14 ) is very stringent, but this is the best we can do for 
this scheme. It is such a stringent assumption that motivated the development of B-j formulation to be 
presented in the next section. 

The proof of this theorem will be given in the following lemmas. 

We define the kernel space 

Xl^ := {77 G X, : b{q,y) = 0, Vy G 

According to the definition, any {u, E, B) G X^^ satisfies (V ■ u,q) = 0, Vg G Qh, and W ■ B = 0. 
Lemma 5. (Boundedness) £i_b(’, •) A bounded on x X^^ and b(-, •) is bounded on Xeh x Y^. 
Proof. We first note that by the energy estimate 

lh“llo, 3 < ||Vrx-|l<i?e|l/||-i, 
and ^ 

||V, xB-||< ||/||_i. 

From Holder inequality and Sobolev imbedding theorem, we see 

|((tl- • V)m,u)| < ||m“||o. 3 ||Vm||||u||o, 6 < ||M“||o. 3 ||M|ll|!lt||l< i?e||/||-l||M||l||r’||l. 

And by Lemma 3 (on X^^), 

|(m X B~,v X B^)\ < ||m X x B^\\ 

<||V,xB-f||n||i||u||i 

and 

|(u X B-,E)\ < ||u X B-||||£;||< IjV;, X B-|||ju||i||£;||< ' ||/||_i|iu||i||£;|l. 

The boundedness of the other terms is obvious. □ 


^ a ^ 


mi sup T—r- II II 

yeYhrK^Xh ll’7llxEl|y||y 

Proof. It is well known that there exists constant ao > 0 such that 

. f -(V ■v,q) 

ml sup ——iTirii” > 0^0 > 0- 


nil sup I, II I 
qeQhvGVh 

This implies for any q G Qh, there exists Vq G Vt, such that 

-(V • Vq,q) > aoljgf , 
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and 

ll«,lli=lkll- 

For the magnetic multiplier, we have V • iFg (div, H) = ^(n)and for any s G there exists 

Cs G -ffg (div, n) such that V • Cg = s and ||Cs||div< (^11511^ where C is a positive constant. 

For any t] = {q, s), choose y = {vg, Cg). Then 

b{v, V) = -(V • Vq, g) + (V • Cs, s) > aolkf+l|sf > min(ao, l)l|y|ly, 

and 

ll^^ll^+IICgllL^ |kf+C2||sf<max(l,C2)|lyf,.. 

This implies the desired result. □ 


The inf-sup condition of a£;(-, •) holds in X^j^\ 


Lemma 7. (inf-sup condition of «£;(•,•)) Assume There exists a positive 

constant (3 such that 

inf sup — > /3 > 0, 

ll^ll-X^E ll^ll-X^E 


and 


inf sup 


aE{^,v) 


ll^llxEUrjIlxE 


> /3 > 0 . 


Proof. We only prove the first inequality, and the second is similar. 
Note that for V • S = 0, there exists i/i G iTg (curl, fl), such that 


V X t/i = —B, 


and 

||t/l||curl<Co||B||, 

where Cq is a positive constant. 

Take v = u, F = E + txl), C = B + RmS~^V x E, where f is a non-negative constant to be 
determined. 

Then 

aE{i,ri) =^||VMf+5|j£; + M X B-f+tS{u x B~,f,) 

itg 

+ tS{E,f,) + tSRf^\\Bf+\\^ X Ef 

From Sobolev imbedding theorem and Poincare’s inequality, we have: 

II«IIo.6<C'i||Vm||. 

Note the following basic estimates: 

M X B^f > ^S\\Ef-S\\u X B~f 

> ^s\\Ef-sc|\\y^, X B-\m\/u\\^ 

> ^S\\E\\^-^R,Cl\\f\\U\\Vu\\f 
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where the last inequality is from the energy estimates. 

X < ^^c^hSR;^^Uf+R^clts\\u X B-f 

< -t-^\\Bf+R^ClClts\\Vh X B-fllVixf 

\ts{E,'ii,)\ < \c^^SR;;,huf+clRraSt\\Ef 
< UsR^^^WBW^+ClRmStWEf. 

Then 

>^||V«f+l5||£;f+ ||V X Ef+\tSR-J\\Bf + ^jVuf 
Taking t = \C^'^R'^, we see 

aE{trj) >^J\Vu\\^+^S\\Er+\\V X Ef+^-C^^SR-^Bf 

By the assumption of R^, we see the last term is non-negative. Hence 

aEitv) > ^J\yu\\^+\s\\Er+\\V X E\\^ + ^Co^SR-^Bf. 

On the other hand, by the choice of ry = {v, F, C), we have 

ll«lli=ll«lii, 

||F||curi< \\E\Ui+\c^^R-^\\xP\Ui< \\E\Ui+^-C^^R-^\\Bl 

\\C\\< ||B||+i?„5-i||£;||curi. 

This implies the desired result. □ 

Under certain conditions, we can show that the Picard iteration converges. We give the detailed proof 
for B-j formulation in §5. With minor modification, the discussion also holds for B-E formulation. 


5. B-J FORMULATION 


As we noted before, the discretization based on the B-E formulation seems not so desirable because 
of the stringent assumption ( 4 . 14 ). In this section, we propose a new finite element scheme whose 
well-posedness will not depend on the assumption such as ( 4 . 14 ). 

We note that it is the variable j that appears in the energy estimate. Therefore it seems natural to use 
B and j as the mixed variables of the electromagnetic part of the MHD system. 

We eliminate E by Ohm’s law and consider the following model: 

— ^Au -f Vp -f (tt • V)tt -f SB X j = f, 

Re 


V X J — V X (m X S) = 0, 


(5-ia) 
(5-lb) 


(5-1C) 

(5-id) 
(5-ie) 
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3 - X B = 0, 

V • M = 0, 

V .B = 0. 

Well-posedness of the continuous formulation has been shown in [15]. The author proved that there 
exists at least one solution u G Hq B G H (curl, fl) n Hq (divO, fl) for the nonlinear system where 
j is eliminated. And now ( 5 . 1 ) is just another form of the system given in [15], if j = V x B is 
considered as an intermediate variable. Therefore we can conclude with the existence of the solution 
of ( 5 . 1 ); for any f G (Tfg there exists at least one solution u G B G i/(curl, fl) n 

iTo(divO, n) and _7 G 

Discretization methods based on B-j formulation actually have already existed in the literature. For 
example, some finite volume methods based on B-j formulation have been developed in [14, 13] where 
the conservation of V • j =0 was considered (but no discussion on the condition V • B = 0), and in 
[16], B and j were used as variables in the simulation of liquid metal breeder blankets. 


5.1. Mixed finite element discretizations of B-j formulation. We now present our new finite element 
discretization of the above system ( 5 . 1 ). 


Problems. Given f G V^. Find {uh, jh,a-h, Bh,ph,rh) G XjhxYh, such that for any {vh,kh,Th,Ch,qh,Sh) G 

^jh X Yfi, 

1 


(5-2a) 
(5-2h) 
f5-2c) 

(5.2d) 

(5-2e) 

f5-2/) 


R, 


-{Vuh,Vvh) + L{uh\Uh,Vh) + S{jh,Vh X Bh) - {ph,^ ■ Vh) = {f,Vh), 

^(V X -^(V X cTh.Ch) + {Th.V ■Ch)=Q. 

-^{cThjTh) - X Bh,Th) = 0 , 

-rtm 

S{jh,kh) - -^{Bh,V X kh) = 0, 

-(V • Uh,qh) = 0, 

(V • Bh, Sh) = 0. 


In the above scheme, an additional variable (Jh is introduced to accommodate for the evaluation of 
the weak curl operator V/jX which is a nonlocal operator. This extra work comes from the nonlinear 
coupling term (V x (it x B), C), because curl operator cannot act onu x B directly. 

Before further discussions, we verify basic properties of the discretization and the energy estimates, 
which are basic and important tools in the design and analysis of numerical methods, especially for 
nonlinear problems. 


Theorem 4. Any solution of Problem 3 satisfies 

(1) Gauss’s law of magnetic field in the strong sense: 

V-Bh=0. 


(2) the Lagrange multiplier rh = 0, hence ( 5 . 2 b) reduces to 

V X {jh - (Th) = 0. 


16 


KAIBO HU AND JINCHAO XU 


(3) the energy estimates: 

-^\\Vuhf+S\\jh\\'^= {f,u), 

and 

+5|b-,f < 

Proof. (1) It is a direct consequence of ( 5 . 2 f), since V • i?g (div; O) = Lg 
(2) Take C?, = V x j/, - V X (T/j. From ( 5 . 2 b) we see 


V X j/, - V X (T/j = 0. 


Hence 

This implies 


[rh, V ■Ch)= 0, yCh G iTo^(div, H). 
T-/! = 0. 


(3) Take Vh = Uh,Ch = Bh,Th =Vh x Bh,kh = jh in ( 5 . 2 a)-( 5 . 2 d). Add them together, we 
have 

■^WyuhW^+SWjhW^+Sijf^Uh X Bh) - X Bh,Vh X Bh) = {f,Uh). 

ilq -t^rn 

Again from ( 5 . 2 d), the last two terms on the left hand side vanish by taking kh = V{uh x Bh). 
This implies the desired result. 

□ 


From ( 5 . 2 c), we see crh — P(M/t x Bh), and from ( 5 . 2 d), we get jh = Rm~^^h x Bh. To prove 
the existence of solution of the nonlinear scheme, we formally eliminate cTh and jh using the above 
identities, to get a system with Uh, Bh, Ph ™d Sh. 

For this purpose, we define 

aj{'il’h-,ih,Vh) ■■= -^C^Uh,\^Vh) + L{wh;uh,Vh) 

ite 

s s 

“f h X Bh,'Oh X Gh) X Gh,^h X Ch) 

g 

+ X BhT^h X Ch). 

Hereafter, -iph, ih, Vh are short for {wh, Gh), {uh,Bh) and {vh, Ch) G Xjh. 

With a slight abuse of notation, we will use the same notation for b{^h, Uh) := —(V • Uh, Qh) + (V • 
Bh, Sh) as in the B-E formulation. But we note that the first variable here l^h = {uh, Bh) G Xjh and 
the first variable ^h = (u, E, B) G Xh in b{^h, Vh) in B-E formulation are in different spaces. 
Eliminating jh and crh. Problem 3 is equivalent to the following form. 

Problem 4. Given 9 = (/, 0) G Xy,find ^h G Xjh, Xh G Yh, such that 

(5-3) aj{^h-,^h,Vh)+ b{vh,Xh) = (9,Vh), W Vh £ Xjh, 

(5-4) b{^h,yh) = 0, yyh&Yh- 

where {9,fih) := {f,Vh). 
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To see the equivalence, we note that if (nit, j/i, cr/j, Vh) € x solves Problem 3, then 

{uh,Bh,Ph,rh) e Xj/.xY/i solves Problem4with the same data and \\{uh,jh,(Th,Bh)\\xy 

Conversely, from a solution (tt/i, of Problem 4, we can reconstruct h x B^^Viuh x Bh),Bh^ph,rh) C 

X Yh which solves Problem 3 with the same data, and 

Wiuh.'^h X X Bh),Bh)\\x^ < ‘^\\{uh,Bh)\\xp 

Such a variational form is closely related to the “curl-formulation”, for example, in [15]. Here curl 
operators are replaced by its discrete version “V/iX”. 

The existence of solution of the nonlinear discrete scheme ( 5 . 2 ) can be stated as 

Theorem 5. There exists at least one solution {uh, Bh,Ph, fh) G Xjh x Y^ solving Problem 4. And 
therefore there exists at least one solution {uh,jh, o’h, Bh,Ph, fh) G Xjh x Yh solving Problem 3. 

It suffices to prove the existence of solution of Problem 4 under the norm 

(5-5) \\{uh,Bh,Ph,rh)n:= ||«,||?+|lB,||^+|b,f+ ||r,f. 

Define the kernel space by 

■= {fth G Xjh ■■ b{rih,yh) = 0, yyh G Yh}. 

Similar to the above discussion, we have 


Lemma 8 . (Boundedness) With the norms given in ( 5 . 5 ), 6 (-, •) is bounded and dj(-; •, •) is bounded in 
From the construction of solutions in Brezzi theory, we note that it is enough to prove the boundedness 

ofdj(-;',-)in^°h- 

Proof. From Cauchy inequality and imbedding theorem, 

{{Uh ■ V)Uh,Vh) < ||M/t||o.3||VM;,|||b?,||o,6< |||| 1 111| M/* || i || i. 

Similarly, 

{{uh ■ V)vh,Uh) < ||M;i||i||||Mii||i|jt;ii|ji. 

On the other hand, from Lemma 3, 

{jh,Vh X Bh) < llV/i X Bh\\\\3h\\c\\vh\\i< \\Bh\\d\\3h\\c\\vh\\i, 

and 

{Uh X Bh,yih X Ch) < ||V;i X Bli||||M;i||i|jC;i||(;< ||Bii||d||M/t|| 1 ||Cii||d. 

The boundedness of other linear terms are obvious. □ 


Here we note again that the estimate of the boundedness of {uh x B/i, Vh X Ch) is a major motivation 
of introducing the modified || -jlc and |j -lb norms, because Uh x Bh may not be in iF^(curl, fl), so curl 
is actually a discrete operator acting on the FIq (div) function Bh. 


Lemma 9. (inf-sup condition of bf, •)) There exists a positive constant a such that 

. r b(fjh,yh) - 

inf sup ——-> a > 0. 

fi^^Xh ll^^llxj lly^ll'5^ 
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Proof. Similar to Lemma 6 , it suffices to prove the following two inf-sup conditions of the pressure and 
magnetic multipliers; there exists constant ao > 0 such that 


inf 

Qh^Qh 


(V • Vh,qh) 


> OiQ > 0 , 


inf sup 

SfeGLg ,,(n) ChGJ/J(div;n) 


(V • Ch,Sh) 


> ao > 0. 


IChWdhhl 

The first inequality is standard for existing Stokes pairs. Now we focus on the second. The proof is a 
three dimensional case of the discussion in Chen et al. [7]. We include the proof here for completeness. 
The major difficulty is that || • is a stronger norm than || • ||div 

It is known that for any Sh € Lg there exists v G i7g such that 

V - v = Sh, 


ll«lli< N^ll- 

Let and 11° be the interpolation in TTg (div, fl) and Lg ^ (we refer to [3] for the definition, and 
[9] for the local bounded cochain projection, which is bounded in i7(curl) and 7T(div)). We denote 
Vh = Then 

v-Vh = v- = n°v • V = n°sh = sh . 

Note that v G iLg (17)°, hence 11°"' is well-defined and bounded. Therefore 

ll«/*ll=l|n°-u||<||n°-||||r;||i< ||n°-||||s^ll- 

Nowit suffices to prove |jV?i x ||s?i||. 

In fact, using inverse inequality and approximation results (see, for example, [ 6 ] and [3]), 

(Vft, X Vh, Vh X Vh) = (Vh X Vh-V X v,Vh X Vh) + (V X V, V/i X Vh) 

= {vh-v,V xVhXVh) + {V xv,VhX Vh) 

< h~'^\\vh - ■w||||V/t X u/i||-f|ju||i||V/i X Vh\\ 

< ||rt||i|!V?, X Vh\\ 

< ||sft||||V;i X Vh\\. 

Therefore 

l|V,.xt;,||<||s^||. 

This proves the desired result. 

□ 


Next we consider to solve the subsystem related to •, •). We introduce the existence theorem for 
nonlinear variational forms, which is given in, for example, [10]. Since we focus on the discrete level 
here, we only given the results for finite dimensional problems. 

Theorem 6. Assume that the dimension of V is finite, and there exists a positive constant a such that 
bounded trilinear form a(-; •, •) on V satisfies 

a{v; V, v) > a\\v\\^, \/vgV. 

Then the problem: given f G V*,find u G V, such that for all v G V, 

a{u-,u,v) = f{v), 
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has at least one solution. 


It is easy to see that 

+IIV;, X B\\\ 

Jte 

From the discrete Poincare inequality (Lemma 4), we have onX°^. Therefore 

the condition in Theorem 6 is satisfied with V = and a(-; •, •) = %(•; •, •). 

Combining Theorem 6 with the boundedness (Lemma 8 ) and the inf-sup condition of b(-, •) (Lemma 
9), we have proved the existence of solution of nonlinear discrete problem (Theorem 5). 


5.2. Picard iteration. In order to solve nonlinear Problem 3, the following Picard iteration can be used: 

Algorithm 2. For n = I,2, i,..., given ^ V,.xiF^(div, F!), / e V*. Find{ul,jl,cTl,Bl,pl,rl) G 


^jh ^ 

Y/i, such that for any [vh, kh 

(Ihi € ^jh 

X Yh, 


( 5 . 60 ) 

^{Vul,Vvn) + L{ul-^ 


sr') 


(5-Qb) 

]^(Vx 



+ {rl,V-Cn) 

( 5 - 6 c) 


-^io-h,Th) - 

s , 

-(n 

Rin 


( 5 - 6 d) 


Sijl.kn, 

)-/- 
-ft 77) 

\Bl,Vxkn) 

( 5 - 6 e) 




-(y-K.Qh) 

(5-6/) 






The following basic properties of Algorithm 2 can be also established similarly. 


Theorem 7. Any solution of Algorithm 2 satisfies 

(1) Gauss’s law of magnetic field in the strong sense: 

y-Bl = Q. 

(2) the Lagrange multiplier = 0, hence ( 5 . 6 b) reduces to 

V X (jX - O = 0. 

(3) the energy estimates: 

^||V<f+5|ljTf= (/,<), 

iLe 

and 

We also recast Algorithm 2 into an abstract form of Brezzi theory for the convenience of analysis. We 
will use to denote , cr^, ) which is the solution of last iteration step (or initial guess). We 

assume uf and are given as known functions. For the initial guess, we assume ||m°||i, ||B°|| and 
lb'°ll= ^ ^h\\ bounded. From the energy estimates, we know \\uf ||i, ||B^ |j and || V/j x B^ || 
are bounded uniformly with the iteration step. 

The variational form can be written as: 
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Problems. Given e Xjh, 0 = {f,l,g,h) e X*^, -ijj = {m,z) G Find {uh, jh,a-h, Bh,ph,rh) G 
Xjh X Yh, such that 

(5-7) aji^h’^h,'nh) + b{r]h,Xh) = {e,r]h), W rjh & Xjh, 

(5-8) b{^h,yh) = {ip,yh), yyh^Yh- 


Now for given •, •) is a bilinear form. Problem 5 corresponds to Algorithm 2 when 

u~ = B~ = B^~^ and I, g,h,m, z = 0. 

We give the main theorem of well-posedness: 


Theorem 8. (Well-posedness of Picard iteration of the B-j formulation) 

There exists unique {uh,jh, CTh, Bh,Ph-, fh) solving Problem 5, and the solution satisfies: 

\\{uh, jh,(Th, Bh)\\x,+\\{Ph,rh)\\Y< C{\\{f ,l,g,h)\\\. + \\{m, z)\\\,). 

where C only depends on the domain and || (/, Z, g, /r)||x*. || (w, z)||y. 

Next we focus on the proof of this theorem. Similar to the nonlinear problem, we first formally 
eliminate the variable jh by x Bh, and formally eliminate (Th to get a system with Uh, Bh and ph, 
Sfi as the variables (Problem 6 below). Boundedness and inf-sup condition of the bilinear form b{-, •) are 
also similar to the nonlinear problem. Finally, we use the coercivity of the bilinear form •, •) on 

Xjf^ to get the well-posedness of the Picard iterations. 

Problem 6. Given G Xjh and e = {f,h) G = {m,z) G Y*,find ^h G Xjh, Xh G Yh, 

such that 

(5-9) aj{if^;^h,'nh)+ bj('nh,Xh) = {9,f]h), W fjh £ Xjh, 

(5-10) bj{ih,yh) = {ip,yh), yyh£Yh. 

where if, Vh) := {f ,Vh) - {l,'P{vh x B")), {h,Ch) ■= {h,Ch) - hX Ch) + {g,V hX Ch)- 


In what follows we use 


to denote the dual norm of FZq (curl, D) (with norm Ij-ljc): 

{l,Fh) 


\l\U:= 


sup 


J’h6ff|f(curl.n) \\-^h\\c 

To see f and h are bounded linear operators, we note the basic estimates: 


(Z,P(n, xB^))<||Z||,.||P(u, xB^)|U 

< |!Z||c*||^t;, X B~\\ 

<\\l\\c*\\vH\\l\\^hXB~\\, 

and 

{l,Vh X Ch) < rnui^h X Ch\\c< ||Z||c*||C;,|U, 

{g,Vh X Ch) < |!g||c*||V„ X Ch\\c< \\g\\c 4 Ch\\d- 

In the following discussion, we will use the Riesz representation lo,go G i7g(curl, D) of l,g G 
FZq (curl, Q)* which are defined by 

{go,Th) := {g,Th), VTh e Hl}{cnrl,n), 


and 


{lo,kh) ■■= {l,kh), Vfe/i G Ffg (curl,fl). 
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Note that ||gol|c= llslU* and ||iollc= M\c*- 

For the relation between Problem 5 and Problem 6, we have; 

Lemma 10. If (u^, Bfi,Ph, fh) solves Problem 6 and 

\\u,,\\l+\\B,\\l+\\p,r+\\r^^^^ 

where C depends on the domain, |l^|l_fr(div)* and\\{m,z)\\Y*, then 

(tr/i, jin Bh, Pin '^h) •— id^h-iBrn ^h ^ Bh d" 5^ Iq, 

^ , Bh, Ph, Th') € ^jh ^ 

solves Problem 5 , and 

(5.11) \\{uh,jh,(Th,Bh)\\\ ^\{ph,rh)\\Y< \\{f ,l,g,h)\W. + \\{m, z)\\'^,. 

^ 3 

On the other hand, if {uh, jh^cTh, Bh,ph,rh) solves Problem 5 , then [uh, BfnPh,fh) solves Problem 

6 . 

Proof In Problem 5 , we take Vh,kh, Ct, Qh, s/t = 0 in (5.7) to see 
(5-12) (Th = P(M?i X Bf^)+ S~^Rmgo, 

and take Vh, Th, Ch,qh,Sh = Ilin (5.7) to note that 

(5-13) jh = X B/i + 

If {uh, BinPh, Th) solves Problem 6, and 

it is easy to see from (5.12) and (5.13) that {uh, hxBh+S-^lo, P(m/ix B~)+S-'^Rmgo, Bh,Ph, rh) 

solves Problem 5 , and 

Rf^WVh X Bn\\l+\\Bn\\l = + I|V • Bnf+{1 + Rf,^)\\Vh x Bnf 

<\\Bh \\1 

||PK X B^)\\< \\un X B^\\< lltt^llillV;, X B^\\. 

We note that || V/i x Bf || is bounded by the right hand side due to the energy estimates. This implies 

\\{uh, R,nV h X Bh + S ^lo,F{uh X Bf^) + S 

<11/1-1 + 11/^11 W + ll(”^’^)ll^* + ll^ll«+ll9ll« 

<\\{f,l,g,h)rx, + \\{m,zWY^. 

3 

On the other hand, solution of Problem 5 also solves Problem 6 by substituting (5.12) and (5.13) into 
( 5 - 7 )- 

□ 

Once the well-posedness of Problem 6 is established, the first part of Lemma 10 will imply existence 
and stability of the original Problem 5 , and the second part will imply the uniqueness. Hence it suffices 
to prove well-posedness of Problem 6 under the norm |j • ^ ((5.5)). 

Similar to the nonlinear case, we have 

Lemma 11. (Boundedness) •, •) is a bounded bilinear form on with respect fo H-Ha ((S-S)) 
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We note that the bound depends on the domain and |jit^||o, 3 , ||V/i x By the energy estimates, 

we know these terms are bounded by known data. 

The boundedness and inf-sup condition of 6 (-, •) are the same as the nonlinear problem. 

Next we show the coercivity of •, •) on 

Lemma 12. There exists a positive constant a such that 

> a{\\uu\\\+\\Bu\\A. G X%. 

Proof. Take Vh = Uh, Ch = Bu, 

aAth-Ah,L) = +^|!V^ X B^f. 

Ite 

From Poincare’s inequality (Lemma 4) and V • Bh = 0 on X^^\ 

||B^||< IIV,, X Bh\\. 

Hence 

||B;,|U< IIV;, X Bh\l 

and there exists a positive constant a which only depends on the domain and R^, Rm, S such that 

□ 

From Lemma 11, Lemma 9 and Lemma 12, we have proved the well-posedness of Problem 6 . And 
from Lemma 10, this shows the well-posedness of Problem 5, and hence Algorithm 2 as a special case. 

5.3. Convergence of Picard iterations. We prove the convergence of Picard iterations for finite ele¬ 
ment B-j formulation under proper conditions. We remark that similar conclusion also holds for B-E 
formulation and the continuous level with minor modification. We omit the subscript “h” in this section. 

Theorem 9. The Picard iteration scheme convergences when 

C!Rl\\f\\.^+ 2 CtRt\\fr_,+ 8 C|RlRU^^^ 

and 

lQRtClRl\\f\W< 1 . 

We remark that the above conditions are satisfied when the data ||/||-i is small relative to R~'^ and 
RfA- This is a common assumption to prove the convergence. 

Proof. By the standard theory of mixed methods, it suffices to consider the convergence in := 
{rjh G Xjh : b{r]h,yh) = 0, Vy/^ S 1^}. 

The n-step equation can be written as 

(5.14) - S{r x B^-\v) = {f,v), 

ILq 

( 5 . 15 ) -(«" X B^-\Wh X C) + X B^,Wh X C) = 0. 

And n — 1-step equation is similarly written as 

( 5 . 16 ) Vt;) - x B^-\v) = (/,«), 

IXq 

X X C) + X X C) = 0. 


( 5 - 17 ) 
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Define the errors 




n—1 


e” :=r-r-\ 


From the equation j" = V/j x B^, we get e" = V/j x e^. 

Subtracting ( 5 .i 6 )-( 5 .i 7 ) from the n-step equation ( 5 .i 4 )-( 5 .i 5 ), we get the error equation; 


1 


1 


1 


- (K-i • V)e-, rt) + - ((er 1 • V)u-\v) - - ■ V)v, el) 

- ^ ((er' • Vr;) + x e^^v) + x = 0, 

X V;, X C) + (e”-i x x C) + x e^, x C) = 0. 

Adding the above two equations and taking ti = e", C = yields 


^((u-^.V)el,el) + ^(iel-^-V)e^ --i 

S' 


) + ;^(Ve;:,Ve:]) + S(er'xj"-\e;:) 


(5.18) + S(e^-' X V„ X C) + — (V„ X e^, Vi, x C) = 0 . 


From the energy estimates, we know 


||VM"|!<Se||/||-l, 


and 

IU"II<(|)’ 

which holds for all n > 0. 

Then we have the estimates for the nonlinear terms; 


^((u^-^-V)el,el) 


- 2 " 

1 

-( 

2 

< 

- 2 1 e 


lo.3||Ve”||||e”||o,6 


II2II ^.n—11 


^ oC'i||Ve"||“||M" ||o ,3 


-illVeJjf, 


((el-^-V)el,n-^) 


< ^IIm 


n—1 1 


o.ellVe^Jlllleri' 


1 0,3 


- 2 1 e 


_i||Ver^||||Ve: 


< j^\\^el-Y+CtRl\\f\\W\^elf, 


|S(er' X r-\el)\< SC^IIV, x eri|||lj"-i||||Ve" 

SSC^S^Ilj—i||||er'lll|VeC|| 


<SC2Rm 


iller'illlVe" 


< ^\\erT+^ReClRl\\f\W\\Wel\\^ 
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and 


Sie^-^ X X e^) < SC2\\Vh x e 


"-1'!I|Vm”- 1||||V;, X e^i 


<5C'2i?ei?^||/||_i||er'lll|e"|| 


< 


,n-l||2 

8K 


e]-Y+8SRl,CiRl\\f\\U\\e-r- 


Combining the above estimates with ( 5 . 18 ), we have 

^ - ic2i?,||/||_i-Cfi?3||^||2_^_4c2^^^^||_^|j2_^ I ||Vgnn2 


+ {j^-SR^^SCiRi\\f\\i, 


|ejf<^||Verif + J-||er'f. 


We define the energy functional to be 




2R ' 

^m. 


Therefore when 


and 




S 


we have 


>8R^^SCiRi\\fr-,. 


gn < -£"-1. 


This implies that (tt", S") converges to some (tt, B) in the norm defined by 

1 ||V«"f + A|iv,xB"||. 

-^e -^rn 

Considering the continuity of the trilinear form, we know the limit can be taken and (u, B) is a solution 
of the nonlinear Problem 3. 

□ 


6. Concluding remarks 

In this paper we considered the mixed finite element discretizations of the stationary MHD system. 
Compared to the time-dependent system, the Gauss’s law of magnetic field is an independent equation 
which cannot be derived from the Faraday’s law. Therefore classical techniques of Lagrange multipliers 
are employed to impose the Gauss’s law. 

Two structure-preserving discretization methods proposed in the previous sections for the stationary 
MHD system preserve both the discrete energy law and most importantly the Gauss’s law V • B = 0, 
but nevertheless they differ significantly in many ways. 

The formulation based on B and E is similar to the time-dependent case studied in [12]. But the 
well-posedness of such a formulation can only be established when the Reynolds number Re is assumed 
to be sufficiently small. To remove such an undesirable constraint, we proposed a new formulation based 
on B and j. Such a formulation was partially motivated by the fact that the energy is given in terms of 
IIjII instead of ||£^||. 

These two formulations look similar. In the finite element discretization of both of the B-E and B-j 
formulations, we have j = E + P(tt x B) (only one variable of E and j is explicitly used in one 
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scheme). This is an equation in i/g (curl, fl). The current density j and the electric field E differ by 
a nonlinear term, which is projected to iJg (curl, ft). But the resulting B-E and B-j formulations are 
different due to the different treatment of the nonlinear term P(m x B) in, for example, the discretization 
of the Lorenz force term. We note that in the B-E formulation, the Lorenz force term {j,v x B) is 
discretized as 

{E u X B,v X B). 

Whereas in the B-j formulation, the corresponding discretization is as follows 

(E + P(m X B), P(u X B)). 

It is easy to see that these two discretizations are indeed different. 

Similar difference can also be found in other places. These differences make the two algorithms 
behave quite differently. In fact, a key point to get well-posedness is the cancellation of the symmetric 
nonlinear coupling terms to get the energy estimates. Under such a restriction, other parts of the schemes 
have to be different according to the different Lorentz force terms. Indeed the energy estimates of 
these two kinds of formulations have already shown the difference. The energy estimates of the B-E 
formulation involve \\E u x BW"^, while the B-j formulation involves |1 j|P= \\E -|-P(m x B)\\‘^. 

As a result of these differences, a careful analysis indicates that the well-posedness of the B-j for¬ 
mulation can be established without any assumption on the size of R^. 
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